<?php

  # HOGFilter project (c) 2009 PV <1minus1@gmail.com>
  # HOGFilter is a free software distributed under the terms of
  # the GNU General Public License http://www.gnu.org/licenses/

  // including required php files
  require_once 'config.php';           // configuration parameters
  require_once 'HOGTree.class.php';    // the Tree class
  require_once 'TreeParser.class.php'; // the TreeParser class

  // reading the species names table from file
  say( "Reading the species names table ... " );
  $bacNamesFile = file( $config['NamesTableFile'] );
  if ( $bacNamesFile === null )
    die( "Can't read names table file '$config[NamesTableFile]'." );
  if ( trim( $bacNamesFile[0] ) !== "abbbrev_hogenom,species_hogenom" )
    die( "Unexpected header line in names table file:\n$bacNamesFile[0]." );
  unset($bacNamesFile[0]);
  $bacNames = array(); // an assoc.array mapping abbreviated species names to their full names
  foreach ( $bacNamesFile as $row )
    if ( ! empty( $row ) )
    {
      list( $abbr, $name ) = explode( ',', trim($row) );
      $abbr = preg_replace( "~^([a-z]+).*$~i", "\\1", $abbr );
      $bacNames[$abbr] = $name;
    }
  unset($bacNamesFile);
  asort($bacNames);
  HOGTree::$bacIDs = array_keys($bacNames);
  say( "done\n" );
  // finished reading the species names table

  // reading trees
  say( "Reading and pre-filtering trees.\n" );
  $treeFilesMask = "$config[TreesDirectory]/*.$config[TreesExtension]";
  $treeFilesList = glob($treeFilesMask);
  if ( empty($treeFilesList) )
    die( "No files found matching mask '$treeFilesMask'." );
  $treeFiles = array();
  foreach ( $treeFilesList as $fileName )
  {
    $treeName = preg_replace( "~\.$config[TreesExtension]$~", "",
                              basename( $fileName ) );
    $treeFiles[$treeName] = $fileName;
  }

  // Beginning of Step 1 of Req.Spec.
  $tn = 0;
  $tc = count($treeFiles);
  $hogTrees = array();
  foreach ( $treeFiles as $treeName => $fileName )
  {
    ++$tn;
    say( "* Input file '$treeName' ($tn/$tc):\n" );
    say( "  Reading... ");
    $newick_tree = file_get_contents($fileName);
    if ( empty($newick_tree) )
      die( "No input." );

    say( "Parsing... " );
    $tParser = new TreeParser($newick_tree);
    $trees = $tParser->getTrees();
    if ( $tParser->errorOccured() )
      die( "Tree parser error: \"" . $tParser->getErrorMessage() . "\".\n" );
    unset($tParser); // don't need a parser any more

    if ( empty($trees) )
      die( "No tree found." );
    $tree = $trees[0];

    say( "Labels... " );
    $hogTree = new HOGTree($tree); // making a tree with parsed labels
    $tree->destroy(); // freeing the memory...
    unset($tree);
    say( "done\n" );

    $skip = false; // if $skip==true -- the tree will be ignored

    $nL = $hogTree->numberOfLeaves();
    $nmin = $config['n'];
    if ( $nL < $nmin )
    {
      say( "  Leaves: $nL < $nmin. " );
      $skip = true;
    }
    else
      say( "  Leaves: $nL. " );

    $nBac = $hogTree->numberOfBacLeaves();
    $nmin = $config['n'];
    if ( $nBac < $nmin )
    {
      say( "Bacterial: $nBac < $nmin. " );
      $skip = true;
    }
    else
      say( "Bacterial: $nBac" );

    if ( !$skip )
      if ( !$hogTree->hasNonBacLeaf() ) // no non-bacterial leaves
      {
        say( " = all. " );
        $skip = true;
      }
      else
        say( ". " );

    if ( !$skip )
    {
      $hogTrees[$treeName] = $hogTree;
      say( "Added to S1." );
    }
    else
      say( "Skipped." );

    say( "\n" );

  }
  say( "Number of trees in S1 - " . count($hogTrees) . ".\n\n" );
  // finished reading trees
  // End of Step 1 of Req.Spec.

  // Beginning of Steps 2 & 3 of Req.Spec.
  say( "Removing trees having p<$config[p]; defining root and removing non-bacterial leaves.\n" );
  foreach ( $hogTrees as $hogName => $hogTree )
  {
    // Step 2 of Req.Spec.
    // filtering by 'p' value
    say( "Tree '$hogName'... " );
    $hogTree->calcD_rec(); // calculating 'd' values in the whole tree
    $dList = $hogTree->getDList();
    rsort($dList);
    list($dm,$dmt) = $dList[0];  // 'd_max' from Req.Spec
    list($dpm)     = $dList[1];  // 'd_predmax' from Req.Spec
    $p = sqrt($dm) - sqrt($dpm); // 'p' from Req.Spec
    say( "p=" . number_format($p,4) . " - " );

    // Step 3 of Req.Spec.
    if ( $p < $config['p'] )
    {
      unset( $hogTrees[$hogName] );
      say( "removed.\n" );
      continue;
    }
    say( "passed. " );
    // finished filtering by 'p' value

    // choosing which half of the tree to remove
    $lc_l  = $dmt->numberOfLeaves();
      // total leaves count in the 'lower' half of the tree
    $blc_l = $dmt->numberOfBacLeaves();
      // bacterial leaves count in the 'lower' half of the tree
    $lc_u  = $hogTree->numberOfLeaves() - $lc_l;
      // total leaves count in the 'upper' half of the tree
    $blc_u = $hogTree->numberOfBacLeaves() - $blc_l;
      // bacterial leaves count in the 'upper' half of the tree
    $chooseUpper = ( $blc_u/$lc_u > $blc_l/$lc_l );
      // choose the upper (== remove the lower) half of the tree ?
    $blc = ( $chooseUpper ? $blc_u : $blc_l );
      // bacterial leaves count in the selected half

    // checking minimal count of bacterial leaves
    if ( $blc < $config['n'] )
    {
      unset($hogTrees[$hogName]);
      say( "Bacterial leaves: $blc < $config[n]. Removed.\n" );
      continue;
    }

    $parent = $dmt->getParent();
    $parent->unlinkChild($dmt);
    $bacTree = ( $chooseUpper ? $hogTree : $dmt );
    $nbTree  = ( $chooseUpper ? $dmt : $hogTree );
    $nbTree->destroy(); // destructing the non-bacterial half to free the memory

    // setting the root to the "breakpoint"
    if ( $chooseUpper )
      $bacTree = $parent;
    $bacTree->becomeRoot(); // not a big deal if the root is already there

    // removing non-bacterial branches
    $bacTree->removeNonBacBranches();

    // putting the tree back to the array (for it's root may have changed)
    $hogTrees[$hogName] = $bacTree;
    say( "Added to S2.\n" );
  }
  say( "Number of trees in S2 - " . count($hogTrees) . ".\n\n" );
  // End of Steps 2 & 3 of Req.Spec.

  // Beginning of Step 4 of Req.Spec.
  say( "Removing edges with low bootstrap values.\n" );
  $resTrees = array();
  foreach ( $hogTrees as $hogName => $hogTree )
  {
    say( "Tree '$hogName'... " );

    $forrest = array();
    array_push( $forrest, $hogTree );
    while ( count($forrest) )
    {
      $curTree = array_pop( $forrest );
      if ( $curTree->numberOfLeaves() < $config['n'] )
        continue; // the tree is too small
      $curTree->removeTransitNodes();
      $bsList = $curTree->getBSList();
      if ( !count($bsList) )
        continue;
      sort($bsList);
      list($mbs,$mbsNode) = $bsList[0];
      if ( $mbs >= $config['b'] )
      {
        $resTrees[] = $curTree;
        continue;
      }
      $tParts = $mbsNode->separate();
      foreach ( $tParts as $tPart )
        array_push( $forrest, $tPart );
    }

    say( "OK\n" );
  }
  // End of Step 4 of Req.Spec.

  // Sending resulting trees to the output file
  $result = "";
  foreach ( $resTrees as $resTree )
    $result .= $resTree->toNewickCompact() . "\n";
  file_put_contents( $config['ResultFileName'], $result );
  // All done.



  // sends the given text to... well stdout, but can send it somewhere else :)
  function say( $text )
  {
    echo $text;
  }

?>
